% Solve model with menu cost or Calvo adjustment friction

totaltic = tic;


%% Parameters 
glob.eta       = eta; 
glob.beta      = betta; 


%% Options
options.cost   = cost; % 1=menu cost, 2=Calvo
options.tolv   = 1e-6; % tolerance on value function (coefficients)
options.minpre = 20;
options.howard = 100;
options.minbreak = 20;
options.minpre   = 100;
options.minbreak = 100;


%% State space parameters;
glob.np        = 100; 
glob.np        = 200; 
glob.nz        = size(Pr_mat,1);
glob.zgrid     = exp(Pr_mat_key)';
glob.Ygrid     = glob.zgrid(:,1);
glob.Pgrid     = glob.zgrid(:,2);
glob.Wgrid     = glob.zgrid(:,3);
glob.Pz        = Pr_mat;
glob.spliorder = 3;


%% State space
glob.pstat     = glob.eta/(glob.eta-1); % static optimal price 
glob.cfstat    = (glob.pstat-1)*glob.pstat^(-glob.eta); % static optimal price 
glob.pmin      = 0.80*glob.pstat; 
glob.pmax      = 1.20*glob.pstat; 
glob.pgrid     = exp(nodeunif(glob.np-2,log(glob.pmin),log(glob.pmax)));
glob.npini     = glob.np;

% Function space and nodes
glob.fspacev   = fundef({'spli',glob.pgrid,0,glob.spliorder},...
                        {'spli',(1:1:glob.nz)',0,1});
glob.fspacev1  = fundef({'spli',glob.pgrid,0,glob.spliorder});
glob.sv        = gridmake(funnode(glob.fspacev));
glob.sv1       = gridmake(funnode(glob.fspacev1));
glob.np        = glob.fspacev.n(1);
glob.pgrid     = glob.sv(1:glob.np,1);
glob.pgrid2    = glob.sv1(1:glob.np,1);
glob.nsv       = length(glob.sv);
glob.nsv1      = length(glob.sv1);

% Monster transition matrices
glob.Pzmonsterv = kron(glob.Pz,speye(glob.np));
glob.Pzminiv    = kron(glob.Pz,speye(glob.np));

% Cash-flow
glob.cf = ( glob.sv(:,1)./glob.Pgrid(glob.sv(:,2)) - glob.Wgrid(glob.sv(:,2)) ) ...
          .* ( glob.sv(:,1)./glob.Pgrid(glob.sv(:,2)) ).^(-glob.eta);

% Initialise guess
try 
    load functions/vold
    cvold  = funfitxy(glob.fspacev,glob.sv,vold);
catch
    vold   = 1./(1-glob.beta)*glob.cf;
    cvold  = funfitxy(glob.fspacev,glob.sv,vold);
end
pAold  = glob.pstat*ones(glob.nz,1);




%% Loop over adjustment costs
if options.cost==1
    xivec = [0 0.000001 0.00001 0.0001:0.0002:0.001 0.002:0.002:0.01 0.05 0.1 0.5 1 5 10]';
    xivec = 1./(xivec/2);
elseif options.cost==2
    xivec = [0.01 0.02 0.05 0.1:0.05:0.3 0.4:0.10:0.9 1]';
end
mean_pA   = zeros(size(xivec));
avg_probA = zeros(size(xivec));

for jj = 1:length(xivec)

    glob.xibar = xivec(jj); 
    
    
    %% 'Bellman' iterations
    
    for i = 1:10000
        % 1. Solve RHS
        % 1.a) adjustment value
        if i<options.minpre || mod(i,options.howard)==0
            obj = @(p) funeval(cvold,glob.fspacev,[p,(1:glob.nz)']);
            pA  = goldenx(obj,glob.pmin*ones(glob.nz,1),glob.pmax*ones(glob.nz,1));
            if glob.pmax-max(pA)<1e-4 || min(pA)-glob.pmin<1e-4; error('bounds bind!'); end
            vA  = obj(pA);
            howard = 0;
        else
            vA  = funeval(cvold,glob.fspacev,[pA,(1:glob.nz)']);
            howard = 1;
        end
        vA    = kron(vA,ones(glob.np,1));
        
        % 1.b) non-adjustment value
        vN  = vold; 
        
        % 1.c) adjustment probability
        xihat = vA-vN;
        if options.cost==1
            xihat = max(0,xihat);
            probA = 1-exp(-glob.xibar*xihat);
            if isinf(glob.xibar)
                probA = ones(size(probA));
            end
        elseif options.cost==2
            probA = glob.xibar*ones(size(xihat));
        end
        
        % 1.d) total value
        if options.cost==1
            v     = - 1/glob.xibar ...
                    + (1/glob.xibar+xihat).*exp(-glob.xibar*xihat) ...
                    + (  probA).*vA ...
                    + (1-probA).*vN;
            if isinf(glob.xibar)
                v = vA;
            end
        elseif options.cost==2
            v     =   (  probA).*vA ...
                    + (1-probA).*vN;
        end
        v     = glob.beta*glob.Pzmonsterv*v;
        v     = glob.cf+v;
        
        % 1.e) new coefficients
        vaux  = reshape(v,[glob.nsv1,glob.nsv/glob.nsv1]);
        cv    = funfitxy(glob.fspacev1,glob.sv1,vaux);
        cv    = cv(:); 
        
        % 2. Norm
        dcv  = norm(cv-cvold)/norm(cvold); cvold = cv;
        dv   = norm(v-vold)/norm(vold); vold = v;
        dpA  = norm(pA-pAold)/norm(pAold); pAold = pA;
        
        % 3. Break?
        if i>=options.minbreak && dcv<options.tolv && dpA<options.tolv;break;end;
    end
    
    if isreal(v)==0; error('Value function contains complex numbers!\n'); end
    
    pAlong = kron(pA,ones(glob.np,1));
    pNlong = glob.sv(:,1);
    
    
    %% Checks
    if dv>10*options.tolv
        error('... VFI algorithm did not converge! \n'); 
    end
    if min(pA)<1.001*glob.pmin
        fprintf('... lowest price close to bottom of grid.\n');
    end
    if max(pA)>0.999*glob.pmax
        fprintf('... highest price close to top of grid.\n');
    end
    
    
    
    %% Save
    save functions/vold vold
    
    
    %% Policy at mean
    
    mean_state = ((glob.nz+1)/2);    %=steady-state mean
    mean_pA(jj) = pA(mean_state);
    
    
    %% Stationary distribution
    % Use policies and transitions P to generate a transition matrix
    % overall state: (kV,xV)
    PzVprime        = glob.Pz(glob.sv(:,2),:);
    fspace_erg      = fundef({'spli',glob.pgrid,0,1});
    pNprime         = funbas(fspace_erg,pNlong);
    pNprime         = dprod(1-probA,pNprime);
    QN              = dprod(PzVprime,pNprime);
    pAprime         = funbas(fspace_erg,pAlong);
    pAprime         = dprod(probA,pAprime);
    QA              = dprod(PzVprime,pAprime);
    Q               = QN+QA;
    test = full(sum(Q,2))-1; if max(abs(test))>sqrt(eps); disp('Something wrong!'); end
    
    % Iterate (fast and more accurate than inverting big horrid matrices)
    
    [mu,d] = eigs(Q',1);
    mu = mu/sum(mu);
    mu(mu<eps) = 0;
    mu = mu/sum(mu);
    % if min(mu)<-sqrt(eps); error('sth wrong'); end
    if min(mu)<-sqrt(eps); 
    
    mu = ones(glob.nsv,1)/(glob.nsv);
    for i = 1:10000
        munew = Q'*mu;
        dmu = norm(munew-mu)/norm(mu);
        mu = munew;
        if dmu<eps
%             fprintf('... i:%d,d: %2.3e\n',i,dmu); 
            break;
        end
    end
    if dmu>eps 
        fprintf('--> Something wrong: stationary distribution did not converge!\n'); 
    end
    mu = mu/sum(mu);
    mu(mu<eps) = 0;
    mu = mu/sum(mu);
    if min(mu)<-sqrt(eps); error('sth wrong'); end
    
    end
    
    test = mu(glob.sv(:,1)==glob.pgrid(1));
    if sum(test)>1e-2; fprintf('--> Something wrong: positive probability in first k-grid-point!\n'); end
    
    % Aggregation
    avg_probA(jj) = mu'*probA;

end


% Average % deviation of price from stst price
avg_dprice  = (mean_pA/glob.pstat-1)*100;







